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Semiconductor  Measurement  Technology: 
TXYZ:  A  Program  For  Semiconductor  IC  Thermal  Analysis 

John  Albers 
Semiconductor  Devices  and  Circuits  Division 
National  Bureau  of  Standards 
Washington,  D.C.  20234 

Abstract 

A  computer  program,  TXYZ,  for  the  thermal  analysis  of  semiconductor  integrated 
circuits  is  presented  and  its  applications  are  discussed.  The  program  makes  use 
of  the  closed  form,  analytic  solution  of  the  steady-state  heat  flow  problem  for  a 
rectangular  three-layer  structure  with  multiple  heat  sources  on  the  top  layer.  The 
temperature  may  be  obtained  for  any  point  or  set  of  points  in  the  structure  and 
is  useful  in  the  determination  of  the  steady-state  thermal  response  of  IC  chips  and 
packages. 

Key  words:  FORTRAN;  Fourier  analysis;  integrated  circuit;  semiconductor;  steady- 
state  heat  flow;  thermal  analysis;  thermal  resistance. 
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EsTTRODUCTION 


Since  the  introduction  of  semiconductor  integrated  circuits,  one  of  the  most  impor- 
tant sources  of  device  failure  is  the  lack  of  temperature  control.  Hence,  an  accurate 
physical  picture  of  the  temperature  distribution  in  the  device  package  under  the 
power  condition  of  actual  operation  is  of  utmost  importance.  The  purpose  of  this 
report  is  to  describe  the  program,  TXYZ,  which  has  been  developed  for  the  thermal 
analysis  of  integrated  circuit  packages.  In  particular,  the  basic  physical  model  and 
the  mathematical  analysis  are  described  and  an  annotated  listing  of  the  program, 
along  with  sample  data,  is  presented. 

The  physical  and  mathematical  model  used  here  is  taken  in  part  from  work  previously 
carried  out  by  Kokkas  [1].  In  order  to  provide  the  reader  with  a  self-contained 
document  for  using  TXYZ,  much  of  the  development  presented  by  Kokkas  has 
been  worked  out  and  specialized  to  the  steady-state  heat  flow  problem.  The  dis- 
cussion which  follows  presents  an  annotated  description  of  Kokkas'  analysis  with 
additional  material  added  where  necessary.  In  particular,  the  equations  have  been 
analyzed  in  detail  so  as  to  investigate  the  convergence  of  the  solutions  used  in 
the  numerical  implementation.  Specifically,  the  form  of  the  solutions  for  small 
and  large  values  of  the  argument  is  shown  to  require  special  consideration  to  avoid 
numerical  overflow  problems. 

This  report  is  naturally  broken  up  into  two  parts.  The  first  deals  with  the  math- 
ematical and  numerical  details  of  the  construction  of  the  program.  The  second 
deals  with  the  program  and  its  use.  For  those  readers  who  are  interested  in  the  use 
of  the  program,  the  section  entitled  "GENERAL  DISCUSSION  OF  THE  TXYZ 
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PROGRAM"  begins  the  portion  of  the  report  where  the  program  and  specific  ex- 
amples are  discussed. 

UNITS 

The  American  semiconductor  industry  has  traditionally  used  mixed  English  and 
metric  units,  but  presently  there  is  a  trend  in  the  direction  of  the  International 
System  (SI)  units.  For  the  purpose  of  conversion,  it  should  be  noted  that  1  mil  = 
0.001  in.  =  25.4  fim  and  that  1  cal/s  =  4.184  W. 
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SOLUTION  OF  STEADY-STATE  HEAT  FLOW  EQUATION: 


SINGLE  RECTANGULAR  LAYER 

Consider  a  material  of  uniform  thermal  conductivity  {ki),  in  the  form  of  a  rectan- 
gular box  of  lateral  dimensions  L^,  Ly  and  thickness  Li .  The  problem  is  to  determine 
the  temperature,  T{x,y,z),  inside  the  material.  It  is  assumed  that  the  temperature 
satisfies  the  steady-state  heat  flow  equation  [2] 

V^T(x,y,z)  =  0.  (1) 

As  this  equation  is  second  order  in  the  three  coordinates,  there  are  six  boundary 
conditions.  Four  of  these  will  be  provided  by  the  lateral  boundary  conditions.  In 
the  present  problem,  all  four  of  the  lateral  boundary  conditions  are  provided  by  the 
assumption  that  there  is  no  heat  flow  out  of  the  lateral  boundaries  of  the  material, 
i.e., 

dx       x=o,L,  dy  j,=o,L„ 

The  remaining  two  boundary  conditions  will  be  provided  by  the  vertical  boundary 
conditions  (in  z).  These  vertical  boundary  conditions  will  not  be  specified  at  the 
present  time  as  the  intent  of  the  present  section  is  to  obtain  a  general  solution  of 
the  one-layer  problem  with  only  the  lateral  boundary  conditions  being  specified. 
As  the  above  equation  is  formulated  in  Cartesian  coordinates,  it  is  convenient  to 
use  Fourier  analysis  techniques  to  solve  the  x  and  y  portion  of  the  equation.  The 
Fourier  transform  with  respect  to  the  variables  x  and  y  is  used,  remembering  that 
the  geometry  is  constrained  to  0,Lx  and  0,  Ly.  This  is  defined  as  [3] 
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ifx,  fy,  z)=         /     T{x,  y,  z)  exp(-2;ri(a;/:,  +  yfy))dxdy,  (3) 
Jo  Jo 


'0  Jo 

where  f^,  fy  are  the  Fourier  transform  variables  which  are  conjugate  to  the  variables 
2,  y.  The  inverse  Fourier  transform  is  defined  as 

•  +  00      /•  +  00 


/•f  oo  r-too 
/       T{fx,  fy,  z)  exp(2;ri(a;A  +  yfy))dhdfy. 
-oo    J  — oo 


(4) 


The  requirement  that  there  is  no  heat  flow  out  of  the  sides  of  the  structure,  i.e., 
dT{x,y,  z)/dx  and  dT(x,y,z)/dy  are  zero  when  x  and  y  are  either  equal  to  zero  or 
to  Lx  or  Ly,  respectively,  leads  to  a  consideration  of  the  expression 

dT(x,y,z)  ^ 
dx 

/+00      p  +  OO  f  N 

J       r(/:c,/j,,^)exp(2;rty/j,)2;r/i<-sin(27r/:ez)  +  icos(2Kfxx)>dfxdfy,  (5) 

and  a  similar  expression  for  dT(x,y,z)/dy.  If  this  is  to  be  zero  at  the  origin,  this 
would  require  that  the  cos  term  be  removed.  Next,  consider  the  resulting  expression 
at  the  other  lateral  boundary,  i.e.,  at  z  =  Z/^  where  it  is  supposed  to  be  zero.  The 
only  way  that  this  could  be  the  case  is  if  the  argument  of  the  sin  function  is  an 
integer  times  tt,  or  that  the  Fourier  transform  variable  is  of  the  form 

=  ^- 

The  same  argument  applies  to  the  y-dependent  portion.  Then,  in  the  Fourier  represen- 
tation, the  temperature  (with  the  above  lateral  boundary  conditions)  may  be  written 
as 

/+00      />  +  oo 
I       T{n,'m,z)cos(n7rx/Lx)cos{in7ry/Ly)dfxdfy.  (7) 
-oo    J  —  oo 
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The  Fourier  transform  equation  is  now  written  as 

r(n,  m,  z)=  /      /     T(x,y,z)cos[nKx/Lx)cos(fmry/Ly)dxdy.  (8) 

As  the  system  is  of  finite  size,  it  is  convenient  to  write  the  integral  in  eq  (7)  as  a 
sum  over  the  Fourier  cos  terms  which  fit  into  the  rectangular  geometry.  Further, 
as  the  cos  function  is  symmetric  around  the  origin,  the  sums  may  be  written  over 
only  the  non- negative  values  of  the  indices.  It  is  important  to  keep  in  mind  that 
the  terms  corresponding  to  m,  n—0  do  not  have  the  factor  of  2  coming  from  the 
symmetry  of  the  cos  function.  In  addition,  the  differentials  may  be  written  as 

df  —  A—  —  ^     ^  _  -'L  —  J_  (g) 

L^r.  Lx  I^X  ^r. 


Then,  the  Fourier  representation  of  the  temperature  may  be  written  as 

rpf        \                  4r(n,  m,  z)  cos{nwx/Lx)  cos{mny/Ly) 
l[x,y,z}=  2^  2^  — — —   .   ,  (lOj 

where  6nn'  is  the  Kronecker  delta  and  is  equal  to  unity  if  n  =  n'  and  zero  otherwise. 
By  substituting  eq  (10)  into  eq  (1)  and  using 


and 
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^  cos{n7rx/Lx)  = —{riTr/Lx)  cos{n7rx/Lx),  (12) 
and  the  same  relation  for  the  y-dependence,  it  is  straightforward  to  show  that 

E°°         4cos(w7r2;/La;)cos(m;rt//L«)  r   ,     ,^        ,      ,^  ^0         \  , 

(13) 


As  the  sum  is  zero  for  arbitrary  values  of  the  variables  x  and  y  (and  the  cos  terms 
are,  in  general,  nonzero),  then  a  necessary  and  sufficient  condition  that  eq  (13)  is 
satisfied  is  that 

|-(n7r/L,)2  -  [rm^lLyf  +  ^^T{n,m,z)  =  0.  (14) 

This  differential  equation  may  be  solved  analytically  using  elementary  methods.  If 
the  variable,  7,  is  defined  as 

then  the  above  equation  may  be  written  in  the  form 

^r(n,  m,  z)  -  i^T(n,  m,  z)  =  0.  (16) 

The  solution  of  this  equation  is 

r(n,  m,z)  =  a  cosh(7>2:)  +  ^  sinh(7>2r),  [17] 

where  the  coefficients  ot  and  ^,  which  may  be  functions  of  7,  are  determined  from  the 
two  -ar-dependent  boundary  conditions.  The  above  equation  is  the  general  solution  for 
the  ^-dependent  Fourier  expansion  coefficients  for  a  single  rectangular  layer.  In  the 
discussion  of  the  problem  of  a  rectangular  y-layer  structure  where  all  of  the  layers 
have  the  same  lateral  dimensions,  the  solution  in  each  of  the  layers  can  be  expressed 
in  the  form  of  the  above  equation  where  the  coefficients  are  to  be  determined  from 
the  two  2r-dependent  boundary  conditions  appropriate  to  each  of  the  layers.  This 
will  be  used  in  the  next  section  where  the  three-layer  problem  will  be  discussed. 
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SOLUTION  OF  STEADY-STATE  HEAT  FLOW  EQUATION: 


THREE-LAYER  RECTANGULAR  STRUCTURE 

The  basic  problem  considered  in  this  section  is  the  calculation  of  the  three-dimensional 
temperature  distribution  in  a  three-layer  structure  which  is  assumed  to  be  of  the 
geometric  form  presented  in  figure  1.  The  three  layers  are  characterized  in  terms 
of  the  thermal  conductivities  and  thicknesses  /c,-,  L,(t  =  1,2,3).  One  particular 
example  of  a  three-layer  structure  is  an  IC  package.  For  this  particular  case,  the 
top  layer  is  the  semiconductor  device,  whereas  the  middle  and  bottom  layers  are 
the  die  attach  and  substrate  layers.  Clearly,  the  thicknesses  and  thermal  conduc- 
tivities of  these  three  layers  are  of  considerable  importance  in  the  dissipation  of  heat 
generated  by  the  power  sources  on  the  surface  of  the  semiconductor  device.  These 
power  sources  are  typically  the  regions  at  or  near  the  surface  of  the  device  where  cur- 
rents are  passed  into  the  device  during  normal  operation.  Consequently,  the  genera- 
tion of  heat  in  the  device  is  one  of  the  unavoidable  side  effects  of  device  operation. 

All  three  layers  of  the  model  are  assumed  to  have  the  same  lateral  dimensions. 
This  assumption  greatly  simplifies  the  analysis  in  terms  of  Fourier  series  which  are 
common  to  all  three  layers.  In  addition,  it  is  assumed  that  there  is  no  heat  flow 
out  of  the  lateral  boundaries  of  the  structure  due  to  either  convection  or  radiation. 
Again,  this  makes  the  analysis  more  manageable. 

The  mathematical  formulation  of  this  problem  is  based  upon  the  following  set  of 
assumptions: 

1)  the  lateral  dimensions  of  the  layers  in  the  structure  are  all  equal  while  the 
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thicknesses  may  be  different, 

2)  there  is  no  heat  loss  from  the  lateral  surfaces  due  to  either  radiation  or  convection 
and  heat  flow  in  the  structure  takes  place  by  conduction, 

3)  there  are  no  heat  losses  due  to  interconnections  to  the  chip  (the  top  layer), 

4)  there  is  no  thermal  contact  resistance  between  the  various  layers  which  are  in 
contact, 

5)  the  heat  sink,  which  is  in  contact  with  the  bottom  layer,  is  ideal  and  has  a 
temperature  equal  to  ambient, 

6)  each  layer  is  of  uniform,  isotropic,  temperature-independent  thermal  conduc- 
tivity, and 

7)  there  is  no  input  power  density  inside  the  structure,  heat  is  generated  only  on 
the  top  surface. 

The  starting  equation  for  the  analysis  is  the  steady-state  homogeneous  heat  flow 
equation  as  discussed  for  the  single-layer  problem, 

V'T(x,y,z)  =  0.  (18) 

In  the  following  discussion,  the  temperature  in  the  first  layer  will  be  referred  to  as 
Ti(x,  y,  z),  while  the  temperature  in  the  second  and  third  layers  will  be  denoted  by 
T2(x,y,z)  and  T2,(x,y,z),  respectively.  Built  into  the  above  assumptions  is  the  fact 
that  the  temperature  and  its  normal  derivative  (proportional  to  the  normal  heat 
flow)  are  continuous  at  the  interfaces  between  the  layers.  The  assumptions  that  heat 
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enters  the  structure  only  on  the  top  surface  through  the  heating  elements  and  that 
the  bottom  of  the  bottom  layer  is  at  the  same  temperature  as  the  heat  sink  are  also 
available.  These  provide  the  six  ^r-dependent  boundary  conditions  on  the  system  of 
equations  generated  by  using  the  solution  of  the  homogeneous  single-layer  problem. 
These  boundary  conditions  may  be  expressed  as  follows.  First,  the  assumption  that 
there  is  no  heat  flow  out  of  the  lateral  boundaries  is  written  as 

dx       x==o,L,  dy  y=o,Ly 

Next,  the  assumption  that  heat  enters  only  where  the  power  is  applied  in  the  first 
layer  is  expressed  as 

OZ  z=0 

where  ki  is  the  thermal  conductivity  of  the  top  layer.  Further,  the  assumption 
that  the  temperature  is  continuous  across  the  interfaces  between  the  layers  may  be 
written  as 

Tt{x,y,z)U^-L,  =  T2(x,y,z)U=.-L,,  (21) 

T2{x,y,z)\^ — (Li+La)  =  Ts{x,y,z)\^ — (Lx+l,)-  (22) 

The  assumption  that  the  heat  flow  is  continuous  across  the  interfaces  between  the 
layers  is  expressed  by  the  conditions  that 

dTi{x,y,z)           _  dT2(x,y,z) 
f^i  ^   —  »   ,  {^^} 

OZ  z  Li  OZ  z=-Li 


dT2{x,y,z)  _  dTi(x,y,z) 

K2   —  K3   • 

dz  z  (L1+L2)  dz  z  (Li+La) 
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Finally,  the  assumption  that  the  temperature  is  continuous  across  the  interface 
between  the  third  layer  and  the  heat  sink  is  written  as 


0 


(25) 


where  all  temperatures  are  measured  relative  to  the  ambient  heat  sink  temperature. 
In  the  above  equations,  k,-  is  the  thermal  conductivity  of  the  i-th  layer  and  Lg  = 
Li  +  L2  +  ^3  •  It  is  also  important  to  note  that  the  origin  of  the  depth  scale  is  at 
the  surface  of  the  top  layer  and  that  all  vertical  distances  are  negative.  For  the 
steady-state  situation,  the  power  density  may  be  written  as 


where  the  function  U[x,  y)  is  the  weighting  function  which  describes  the  geometry 
and  uniformity  (or  nonuniformity)  of  the  heat-generating  components,  and  Pq  is 
the  steady-state  power  density  per  unit  area.  As  discussed  previously,  the  ambient 
temperature  trivially  satisfies  the  heat  flow  equation.  Hence,  it  is  convenient  to 
consider  all  temperatures  relative  to  ambient.  Therefore,  when  the  power  density  is 
set  equal  to  zero,  the  temperature  will  be  zero  as  it  is  relative  to  ambient.  As  the 
three  layers  are  of  the  same  lateral  dimensions  and  there  is  no  heat  flow  out  of  the 
lateral  boundaries,  the  temperature  in  each  of  the  three  layers  may  be  written  in 
the  form  of  eq  (10).  Further,  the  ^-dependent  Fourier  expansion  coefficients  will  be 
of  the  form  of  eq  (17)  and  may  be  written  as 


(26) 


ri(n,  m,  z)  =  ai  cosh(7xr)  +  sinh(72?), 


(27) 


T2{n,  m,  z)  —  a2  cosh(7>2r)  +  ^2  sinh(7z). 


(28) 
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Ts{n,  m,  z)  —  0:3  cosh(72r)  +  sinli(72r). 


(29) 


There  are  six  unknowns  involved  in  these  solutions  and  they  are  specifically  the 
set  of  expansion  coefficients  a{,^i{i  =  1,2,3).  These  coefficients  may  be  explicitly 
evaluated  by  means  of  the  boundary  conditions  on  the  temperature  and  its  deriva- 
tive evaluated  at  the  interfaces.  The  Fourier  coefficients  may  be  used  directly  in 
the  above  boundary  conditions  as  the  temperature  is  a  sum  over  the  Fourier  expan- 
sion coefficients.  By  substituting  the  above  equations  into  the  appropriate  boun- 
dary condition  equations,  it  is  possible  to  obtain  a  system  of  six  equations  in  six 
unknowns.  For  the  present  case,  this  system  reduces  to  one  with  four  equations 
in  four  unknowns.  The  standard  method  for  solving  this  system  is  to  use  Cramer's 
rule  [4]  with  the  Laplace  method  [4]  for  the  evaluation  of  the  determinants  involved. 
However,  instead  of  using  this  method  explicitly  here,  the  Fourier  coefficients  for 
each  of  the  layers  will  be  presented  and  will  be  shown  to  satisfy  the  heat  flow  equa- 
tion and  the  appropriate  boundary  conditions.  In  particular,  the  Fourier  coefficients 
in  the  three  layers  are  given  by  eqs  (15-17)  in  Kokkas'  paper.  Specializing  these  to 
the  steady-state  situation,  these  are 


72 (n,  m,  z)=a\d  cosh(7(Li  -i- L2  +  z)) E sinh(7(Li  -h  L2 z))  [,  (31) 


(n,  m,  z)  =  AIB  cosh(7(Li  -\-  z))  -\-  C  sinh(7(Li  +  z)) 


(30) 


(32) 


where 


(33) 
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B  =  D  cosh(7Z/2)  +  E  smh(7L2), 


(34) 


C  =  sinh(7L2)  4-  £;cosh(7L2)|, 


(35) 


D  =  siiih(7L3),  (36) 


E  =      coshir^Ls),  (37) 

K>2 


U{n,m)=  I      I     U{x,y)cos{n7rx/Lx)cos[m7ry/Ly)dxdy,  (38) 
^0  Jo 


and 

0  Jo 

is  the  double  Fourier  cos  transform  of  the  power  density  uniformity  function.  Now, 
it  will  be  shown  that  the  above  are  the  solutions  of  the  equation  (see  eq  (16)) 

Ti(n,  m,  z)  -  7^  n(n,  m,  z)  =  0,  (39) 

where  the  subscript  i  takes  on  the  values  of  1,2,3.  This  may  be  easily  shown  to  be 
the  case  as 


and 


Q2 

cosh(7(L  +  z))  =  7^  cosh(7(L  +  z)),  (40) 


q2 

sinh(7(L  +  z])  =  7^  sinh(7(L  +  z)),  (41) 


dz^ 

where  L  is  a  constant  and  is  equal  to  Li,  Li  +  L2,  or  Lg  in  eqs  (30-32).  Hence,  eqs 
(30-32)  satisfy  the  2r-dependent  differential  equation.  Next,  it  will  be  shown  that 
these  Fourier  coefficients  satisfy  the  appropriate  boundary  conditions.  The  first  of 
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these  is  that 


«i  '  ^        =  P(n,  m)  =  U(n,  m)Po .  42 

az  z=o 

Using  the  Fourier  coefficient  given  by  eq  (30)  for  the  top  layer,  this  may  be  evaluated 
as 

^7j_(n^_r7V^I       _  ^^^^]  Q  _|_  (7  cosh(7Zyi)| 

oz         z=0  {  J 

/ci7      i^i?sinh(7Li)  +  C/cosh(7Li)  J  J 


=  U{n,m)Po.  (43) 

Hence,  the  top  layer  boundary  condition  is  satisfied  by  the  ri(n,  m,  z).  Next,  consider 
the  bottom  layer  boundary  condition,  i.e., 

Tz[n,m,z% — =  0.  (44) 

Making  use  of  eq  (32),  this  may  be  readily  evaluated  as 

Tz(n,  m,  z%=-L,  =  Asinh(7(L^  -  Lg))  =  0.  (45) 

Hence,  the  last  boundary  condition  is  satisfied.  The  final  set  of  boundary  conditions 
to  be  verified  are  the  ones  which  pertain  to  the  interface  boundary  conditions.  These 
are 

Ti(n,m,z% — Lr=  T2(n,m,z%=.-Li,  (46) 
T2(n,  m,  z)\^ — (L1+L2)  =  Ts{n,  m,  z)\^ — (L1+L2),  (47) 
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dTi{n,m,  z)  dr2(n,m,z) 

dz  L=-L.  =   dz  U-L/ 

dT2{n,m,  3r)  drs(n,m,  z)  ^ 

dz  z=-{Li+L2)  dz  z^-{Li+L2) 

In  order  to  verify  these  equations,  it  is  simplest  to  calculate  the  right-  and  left-hand 
sides  of  the  equations  and  then  compare  them  directly.  The  left-hand  side  of  the 
first  temperature  continuity  equation  may  be  evaluated  as 

Ti{n,m,z)\^ — =  Tlj^ cosh(7(Li  -  Li))  +  Csinh(7(Li  - 

=  AB  =  A|z)cosh(7L2)  +  J^sinh(7L2)|-  (50) 
The  right-hand  side  of  the  equation  may  be  evaluated  as 

T2{n,m,z)\^ — Li  =         cosh(7(Li  +  L2  -  ^i))  +  ^sinh(7(Li  H-Z^  -Li))| 

=  aI^D  cosh(7L2)  +  E  sinh(7L2)|.  (51) 

Hence,  the  first  of  the  temperature  continuity  equations  satisfies  the  boundary 
condition.  Next,  consider  the  second  temperature  continuity  equation. 

r2(n,  m,  z)\^ — (L1+L2)  =  rs{n,  m,  z)\^ — (L1+L2),  (52) 
The  left-hand  side  of  the  equation  may  be  evaluated  as 

T2(n,m,z)\^  (L1+L2)  =  A|Dcosh(7(Li+L2-Li-L2))+£;sinh(7(Li+L2-^i-^))| 


=  AD  =  Asinh(7Z/3). 
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(53) 


The  right-hand  side  is 

rs{n,  m,  z% — {L^^■L^)  =  Asinh(7(L^  -  Li  -  L2))  ==  Asinh(7L3).  (54) 

Hence,  the  second  temperature  continuity  equation  is  satisifed.  The  next  equation 
to  be  verified  is  the  first  heat  flow  continuity  equation,  i.e., 

dTi{n,m,  z)  dT2(n,m,z) 

 dz  1       r   dz  1       t'  ^^^^ 

OZ  z=—Li  OZ  z=—Li 

Evaluating  the  left-hand  side  leads  to 

=  =  AC27^|^sinh(7L2)  +  £;cosh(7L2)|-  (56) 

The  right-hand  side  may  be  evaluated  as 


iji)  sinh(7L2)  +  E  cosh(7L2)|. 


^^dT2(nm,z)^  =  /C27A^i)sinh(7L2)  +  ^cosh(7L2)^  (57) 

uZ  z=—Li 


The  final  equation  to  be  evaluated  is  that  for  the  heat  flow  continuity  between  the 
second  and  third  layers,  i.e., 

dz  z  (L1+L2)  dz  z=-{Li-j-L2) 

The  left-hand  side  of  this  equation  may  be  evaluated  as 

dT2(n,m,z) 


K,2 


dz  z=-{Li+L2) 


=  /c27a|d  sinh(7(Li  -\-  L2  -  Li  -  L2)) E  cosh(7(Li  +  L2  -  Li  -  L2))| 
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(59) 


The  right-hand  side  is 


dn{n,m,z 


=  Ki^Acosh['^(Lz  —  Li  —  L2)  =  /C37j4.cosh(7Z/3).  (60) 


dz 


z=-{Li+L2) 


Now  that  it  has  been  shown  that  the  Fourier  coefficients  satisfy  the  steady-state 
heat  flow  problem  and  the  appropriate  boundary  conditions,  there  are  several  points 
to  be  considered  before  getting  into  the  body  of  the  program.  These  include:  (1) 
evaluation  of  the  function  U(n,m)  for  a  uniform  power  source  of  given  size  and 
(2)  simplification  of  the  Fourier  coefficients  for  subsequent  numerical  analysis.  The 
latter  point  is  necessary  as  the  limits  of  7  very  small  and  7  very  large  may  give  rise 
to  overflow  or  underflow  problems  when  the  program  is  constructed. 
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FORM  OF  THE  FUNCTION  U(n,m) 

The  first  thing  to  be  considered  is  the  form  of  the  function  U(n,  m)  for  an  arbitrary 
number  of  heat  sources.  This  function  is  defined  as 


The  analysis  can  most  easily  be  accomplished  in  terms  of  a  single  uniform  heat 
source.  The  case  of  an  arbitrary  number  of  heat  sources  follows  by  summing  the 
results  of  each  heat  source.  Further,  if  any  of  the  heat  sources  are  nonuniform, 
their  effects  can  be  constructed  by  suitably  overlapping  a  number  of  uniform  heat 
sources.  In  the  coordinate  system  being  used,  consider  a  single  heat  source  denoted 
by  the  index  i  with  a  corner  at  the  location  [xi,  t/,)  and  lengths  along  the  x-  and 
y-directions  given  by  (lxi,lyi).  Over  the  area  of  the  heat  source,  U(x,y)  is  assumed 
to  be  uniform  and  equal  to  unity.  Away  from  the  area  of  the  heat  source,  U(x,y) 
is  assumed  to  be  equal  to  zero.  Then,  ?7(a;,  y)  may  be  viewed  as  being  a  unit  step 
function  over  the  surface  of  the  power  source.  Consequently,  the  contribution  from 
this  single  heat  source  may  be  written  as 


^0  Jo 


(62) 
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The  integrals  can  be  simply  evaluated  to  give  the  result  that 


Ui{n,  m)  = 


(n7r)(m;r) 


X 


(63) 


Similar  expressions  may  be  "written  for  each  of  the  heat  sources  and  then  summed 
to  give  the  cumulative  heat  source  effect.  In  addition,  as  indicated  previously,  if 
there  are  any  nonuniform  heat  sources,  their  effect  can  be  constructed  by  means  of 
overlapping  uniform  heat  sources.  Before  turning  to  the  small  7  and  large  7  behavior 
of  the  Fourier  coefficients,  it  is  important  to  consider  the  behavior  of  the  function 
U{n,m)  for  either  n  =  0,  m  =  0,  or  both.  This  is  important  in  the  numerical 
implementation  of  the  solutions  as  the  program  will  have  to  calculate  the  double 
Fourier  cos  transform  over  the  range  of  n,  m  required  by  the  sum  in  eq  (10).  Once 
the  value  of  n  or  m  is  zero,  there  will  be  problems  with  most  machines  as  far  as 
evaluating  the  seeming  divergence.  This  can  be  circumvented  by  investigating  the 
behavior  of  the  function  for  n  =  0.  This  may  be  readily  carried  out  by  first  noting 
that  the  function  U(n,  m)  (for  the  particular  form  of  U(x,  y))  is  a  product  of  two 
terms.  This  may  be  simply  written  as  U{n,m)  =  U{n)U{m).  Then,  consider  the  n 
(from  the  x  integration)  contribution  to  the  function  which  is  given  by 


There  is  an  apparent  divergence  or  infinity  if  n  is  simply  set  equal  to  zero.  This  is 
the  way  in  which  a  computer  would  look  at  the  expression.  However,  this  infinity  is 
not  real  as  can  be  seen  by  using  the  expansion  of  the  sin  function  for  small  values 


(64) 
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of  the  argument.  In  particular 


sin(a;)  =  z  -  —  +  •  • 


(65) 


Making  use  of  this  expression  for  the  sin  function,  it  is  straightforward  to  show  that 


The  same  conclusion  holds  for  the  y-dependent  portion,  i.e.,  U{m).  This  must  be 
specially  coded  to  bypass  any  overflow  problem.  The  specific  coding  of  the  heat 
source  may  be  found  in  the  program  listing  in  the  function  UZERO. 

In  general,  the  function  U{n,  m)  is  oscillatory  and  does  not  approach  zero  sufficiently 
fast  for  large  values  of  n  or  m.  In  particular,  eq  (63)  shows  that  Ui{n,m)  0 
Uke  1/nm  as  m,n  ->  oo.  In  addition,  the  cos  terms  in  eq  (10),  i.e.,  cos[nKx/Lx) 
cos{m7ry / Ly),  do  not  approach  a  definite  limit  for  large  values  of  the  argument. 
Consequently,  the  product  U{n^7n)cos(nKx/Lx)cos{mwy/Ly)  tends  to  zero  slowly. 
An  example  of  the  behavior  of  the  z-dependent  portion  of  this  product  function, 
i.e.,  U{n)  cos{n7rx / Lx)  is  presented  in  figure  2.  The  heat  source  used  for  this  figure 
is  a  1  mil  by  1  mil  heat  source  centered  on  the  surface  of  a  200  mil  by  200  mil 
structure.  The  product  function  is  evaluated  at  the  midpoint  of  the  heat  source. 
Because  of  the  symmetry  used,  the  function  shows  damped  oscillatory  behavior.  It 
is  clear  from  the  figure  that  the  product  function  does  not  fall  off  to  zero  sufficiently 
fast  and  hence  does  not  provide  for  rapid  convergence  of  the  calculated  temperature 
using  the  Fourier  representation. 


lim  U{n)  =  lim  — 
re— 0  re— 0  nw 


(66) 
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BEHAVIOR  OF  FOURIER  COEFFICIENTS: 


SMALL  VALUES  OF  THE  ARGUMENT 

As  lias  been  seen  in  the  treatment  of  the  function  U{n,m),  care  must  be  taken  for 
the  case  where  both  n  and  m  are  equal  to  zero  (or  7  =  0).  The  same  considerations 
must  be  carried  out  for  the  Fourier  coefficients  in  the  three  layers.  In  the  solutions 
in  the  three  layers,  the  summation  indices  n,  m  appear.  The  summation  over  these 
variables  is  of  influence  in  the  variable  7  according  to  eq  (15).  Also,  the  Fourier 
coefficients  contain  the  hyperbolic  functions  which  depend  upon  7.  For  large  values 
of  7,  the  sinh  and  cosh  functions  grow  exponentially.  This  can  present  special 
numerical  problems  when  the  summation  variables  approach  the  upper  limits  which 
may  be  required  for  the  case  of  very  small  heat  sources.  Hence,  special  care  must  be 
taken  to  study  the  behavior  of  the  Fourier  coefficients  for  small  7  and  large  7  so  as 
to  remove  any  potential  numerical  overflow  problems.  Once  this  is  properly  taken 
care  of,  the  Fourier  coefficients  and  the  solutions  will  be  numerically  well  behaved. 

First,  consider  the  small  7  behavior  of  the  Fourier  coefficients.  This  is  done  by 
considering  the  small  7  behavior  of  the  eqs  (30-37)  and  the  small  argument  behavior 
of  the  hyperbolic  functions.  In  the  following  discussion  as  well  as  the  discussion  of 
the  large  7  behavior,  the  term  U(n,m)Po/ k.i  will  be  removed  for  convenience.  This 
term  will  henceforth  be  included  explicitly  in  the  sum  in  eq  (10)  as  the  Fourier 
coefficients  for  all  three  layers  contain  this  as  a  common  factor  through  A  (see  eqs 
(30-33)).  Then,  for  small  7, 

Ef^'^,  (67) 

K2 
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Df^^Ls,  (68) 


(69) 


Bf^qLs-h  —1L2,  (70) 

and,  remembering  that  the  factor  U{n,m)Po/i<ii  has  been  included  explicitly  in  eq 
(10), 

A^-'^.  (71) 

Making  use  of  these  expressions  and  the  small  argument  behavior  of  the  hyperbolic 
functions,  it  is  straightforward  to  investigate  the  small  7  behavior  of  the  solutions. 
In  particular, 

ri(n,m,z)  ^  -  —  Ll^  +  —1L2  +  —(7% ^3  +  —\i(Li-\-z)\ 

7  /C3  [  K2  K,i(  K2}  J 


— <  L3  +  — L2  +  — (m  +  ^)  r 

t  K2  J 


Then, 


lim  Ti(n,  m,  z)  =  (Li  +  2r)  +  — L2  +  —  L3.  (72) 

7-*0  /C2  /C3 


Next, 


T2(n,m,z)      i  — I7L3  +  —7(^1  +1^  +^^)1, 

7  «3  I  K2  J 


or 


lim  r2(n,  m,  ^r)  =  — L3  +  ~{Li  +  L2  +  -sr).  (73) 

7-*0  ^3  K2 
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And  finally, 


Ts{n,m,z)      -  —  \-f(Li  +  L2  +     +  z)\, 


or 


lim  Ts{n,  m,z)  =  — (Li  +  L2  +  -^'3  +  z).  (74) 

These  special  forms  of  the  Fourier  coefficients  (in  the  limit  as  7  ->  0)  are  necessary 
in  the  code  to  bypass  overflow  problems  for  small  values  of  the  argument. 
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BEHAVIOR  OF  FOURIER  COEFFICIENTS: 
LARGE  VALUES  OF  THE  ARGUMENT 

As  the  Fourier  coefficients  have  been  investigated  for  small  7  and  have  been  shown  to 
be  well  behaved  when  properly  written,  what  remains  is  to  write  these  coefficients 
in  a  form  which  is  amenable  for  investigating  their  large  7  behavior.  As  noted 
before,  the  hyperbolic  functions,  sinh  and  cosh,  grow  exponentially  for  large  values 
of  the  argument.  On  the  other  hand,  the  hyperbolic  tanh  approaches  unity  for  large 
values  of  the  argument.  With  this  in  mind,  let  us  investigate  the  form  of  the  Fourier 
coefficients,  written  as  much  as  possible  in  terms  of  the  tanh,  which  takes  care  of 
this  potential  numerical  difficulty.  To  this  end,  it  is  convenient  to  introduce  the 
shorthand  notation  for  the  hyperbolic  functions,  c(x)  =  cosh(a:),  8{x)  =  sinh(a;), 
and  t{x)  =  tanh(a;).  Making  use  of  this  shorthand  notation,  the  Fourier  coefficients 
may  be  written  as 


(75) 


(n,  m,  z)  =  A\  Dc{^(Li  +  L2  +  z))  +  Es^^iL^^  +  L2  +  -sr))  [  (76) 


Ts{n,m,z 


)  =  A8(q(L^  +  0). 


(77) 


where. 


(78) 


B  =  I>c(7L2)  +  Es(^L2), 


(79) 
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C  ^  —  D8(^L2)  +  Ec{^L2)  , 


(80) 


(81) 


(82) 


As  in  the  investigation  of  the  small  7  behavior  of  the  Fourier  coefficients,  the  factor 
U[n,  in)Po/Ki  has  been  deleted  from  eq  (78)  for  convenience.  This  factor  may  simply 
be  included  in  the  Fourier  representation  of  the  temperature,  eq  (10),  as  it  is  common 
to  all  three  layers.  Now,  the  above  equations  (eqs  (75-82))  will  be  rewritten  by 
making  use  of  the  definition  of  the  hyperbolic  tanh,  i.e.,  t(x)  =  3{x)/c(x).  First, 
consider  the  coefficient  C. 


C  =  — 


c 


(83) 


Next,  the  coefficient  B 


may  be  written  as 


B  =  Dc{'^L2)  +  ^5(7^2) 


B  =  8{r^Ls)c(^L2)  +  ^c(7L3)«(7L2) 


(84) 
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Then, 

B5(7Li)  +  Cc(7Li)  = 
=  c(^LMnL2)c[iU){t[^L^)t[^^^^  (85) 

^  AC2  K.I  Kl  ) 

Then,  the  coefficient  A  may  be  written  as 


(86) 

It  is  convenient  to  define  the  function  0(7)  as 

n(7)  =  -,  (87) 

{«(7i,)<(7Li)  +  =f  t(iZ„)%I^)  +  "^ti^LM^L,)  +  ^1 
which  is  well  behaved  for  all  values  of  7.  Then,  the  coefficient  A  may  be  written  as 

_  "(7)  (88) 

7c(7L3)c(7L2)c(7A)' 
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By  making  use  of  the  above  procedure,  it  is  relatively  straightforward  to  show  that 
the  Fourier  coefficients  as  given  by  eqs  (75-77)  may  be  written  as 


Ti  (n,  m,  z)  = 


n(7)<7(^i  +  z)) 


72  (n,  m,z)  = 


and 


73  (n,  m,  z)  = 


7c(7Li)c(7L2)c(7L3) 


(91) 


In  eqs  (89-91),  the  function  0(7)  and  the  terms  inside  the  curly  brackets  are  well 
behaved  for  all  values  of  the  variable,  7.  The  sinh  and  cosh  terms  which  remain 
may  still  give  rise  to  numerical  overflow  problems  for  large  values  of  the  argument. 
However,  as  both  of  these  functions  grow  exponentially  for  large  values  of  the 
argument  and  they  appear  in  both  the  numerator  and  the  denominator  of  the 
Fourier  coefficients,  there  will  be  cancellation.  This  cancellation  for  large  values  of 
the  argument  will  not  be  worked  out  in  detail  here  but  is  contained  in  the  FORTRAN 
listing  of  the  program  in  the  function  FUNZ.  These  Fourier  coefficients  are  used  in 
the  equation 


for  i  =  1,  2,3  to  obtain  the  solutions  in  each  of  the  three  layers.  In  the  above,  the 
term  PoU{n,m)/Ki  has  been  written  out  explicitly  and  is  no  longer  contained  in 


00  00 


4U{n,  m)r,(n,  m,  z)  cos{nwx/Lx)  cos{m7ry / Ly) 

i^nO  +  l)((^mO  +  l)L:^LyK,i 


Ti{x,y,z)  =  Po  J2  E 


(92) 


«==0  m=0 
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the  Fourier  coefficients.  It  is  important  in  obtaining  the  solution  in  a;,  y,z  to  use  the 
appropriate  layer  equation.  This  is  automatically  taken  care  of  in  the  program  as 
the  depth  z  is  compared  with  the  various  thicknesses  and  the  corresponding  layer 
Fourier  coefficient  is  used.  The  user  does  not  have  to  specify  which  layer  is  to  be 
used. 

An  interesting  exercise  left  to  the  reader  is  to  show  when  the  three  thermal  conduc- 
tivities are  equal  that  the  one-layer  solution  is  obtained.  Also,  another  exercise  is 
to  show  that  the  one-layer  solution  is  obtained  when  the  thicknesses  of  the  second 
and  third  layers  are  set  equal  to  zero. 

In  the  subsequent  discussion,  it  is  important  to  keep  in  mind  that  the  Fourier 
coefficients  are  functions  of  the  variables  n  and  m.  As  discussed  in  the  previous 
sections,  the  power  density  function  U{n,  m)  tends  to  zero  very  slowly  for  large 
values  of  the  argument.  Also,  the  cos  terms  in  eq  (92)  do  not  tend  to  any  limit  as  the 
arguments  approach  infinity.  It  is  the  Fourier  coefficients  which  are  responsible  for 
the  convergence  of  the  sum.  Consequently,  some  discussion  of  the  n  and  m  behavior 
of  the  Fourier  coefficients  is  warranted.  To  show  how  the  Fourier  coefficients  behave 
as  a  function  of  n  and  m  (or  7),  these  functions  were  studied  in  some  detail.  In  figure 
3,  a  three-dimensional  plot  of  the  n  and  m  dependence  of  the  top  layer  Fourier 
coefficient  is  shown  at  z  =  0  for  a  three-layer  structure.  In  the  figure,  it  can  be  seen 
that  the  Fourier  coefficients  are  peaked  at  the  origin  of  the  n,  m  plane  and  approach 
zero  as  n  and  m  become  large.  In  addition,  figure  4  contains  the  same  type  of  plot 
for  the  situation  where  all  three  layers  have  the  same  thermal  conductivity,  i.e.,  for 
a  thick  one-layer  structure. 
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The  dependence  on  n  of  the  three  Fourier  coefficients  at  the  top  of  each  layer  is 
contained  in  figure  5  for  the  structure  used  to  generate  figure  3.  The  curve  marked 
by  A  shows  the  behavior  on  the  top  surface  of  the  structure,  i.  e.,  at  z  =  0.  The 
curves  denoted  by  B  and  C  present  the  results  for  the  top  surfaces  of  the  second 
and  third  layers,  respectively.  Figure  6  is  a  log-log  plot  form  of  figure  5a  and  is 
especially  important  as  it  shows  the  slow  convergence  of  the  top  surface  Fourier 
coefficient  which  will  be  discussed  in  detail  in  the  section  on  the  effects  of  the  upper 
limit  on  the  calculated  temperature. 
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SPECIAL  CASE  OF  POWER  SOURCE  COVERING  TOP  SURFACE 


As  a  special  case  of  eq  (92),  consider  the  situation  of  a  single  power  source  completely 
covering  the  top  surface;  i.e.,  there  is  a  single  heat  source  with  lateral  dimensions 
equal  to  that  of  the  three-layer  structure.  In  this  particular  case,  it  will  be  shown 
that  the  above  equation  reduces  to  the  familiar  thermal  resistance  equation.  The 
easiest  way  to  proceed  with  the  analysis  is  to  consider  the  specific  form  of  the 
function  U{n,m).  From  eq  (63), 

Ui  (n,  m)  = 


(93) 

For  the  particular  case  of  uniform  surface  coverage,  xi  =  yi  =  0,  Ixi  =  Lx  and 
lyi  =  Ly.  Upon  substituting  these  values  into  the  equation,  the  function  reduces  to 

Ui{n,  m)  =  sin(m7r)|.  (94) 

This  is  zero  when  the  indices  are  nonzero.  For  the  case  where  both  of  the  indices 
are  zero,  the  use  of  the  expansion  of  the  sin  function  gives  rise  to  the  result  that 


Ui(n,m)  =  LxLy6nQ6mQ-  (95) 

If  this  form  of  the  U[n,  m)  function  is  substituted  into  the  equation  for  the  tem- 
perature in  each  of  the  three  layers  (eq  (92))  and  the  form  of  the  Fourier  expansion 
coefficients  as  7  — ►  0  (eqs  (72-74))  is  used,  it  is  readily  shown  that  the  temperatures 
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in  the  three  layers  may  be  written  as 

r.(.,!,,^)  =  Po{^  +  ^  +  ^}  (95) 

Ki  K2       K,^  ) 

Ux,y,.)^Po{'^+^  +  '  +  ^]  (96) 

Ts(x,y,z)  =  Fo|  y  (97) 

As  Pq  is  the  power  density  per  unit  area,  these  equations  give  rise  to  the  usual 
results  of  the  one-dimensional  calculations  of  the  thermal  resistance. 
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EFFECT  OF  UPPER  SUMMATION  LIMITS  ON  TEMPERATURE 

Previously,  the  dependence  of  the  Fourier  coefficients  on  n  and  m  has  been  discussed. 
Figures  3  to  6  contain  typical  results  of  this  behavior.  The  purpose  of  the  present 
section  is  to  show  how  this  behavior  is  mirrored  in  the  calculated  temperature.  In 
particular,  eqs  (89-91)  will  be  used  in  eq  (92)  with  a  variable  upper  summation  Hmit 
in  the  latter  equation.  To  simplify  the  analysis,  a  single  heat  source  with  y,-  =  0 
and  Ij/i  =  Ly  will  be  used.  The  width  of  this  stripe  heat  source  will  be  Ixi  —  1 
mil.  Finally,  the  lateral  dimensions  of  the  three-layer  structure  will  be  taken  as 
Lx  =  Ly  =  200  mil.  This  particular  choice  will  require  only  the  m  =  0  term  in 
the  sum  to  be  retained  while  needing  a  large  value  of  n  terms  in  order  to  obtain 
convergence  of  the  sum.  As  there  is  complete  coverage  along  the  y-direction,  the 
7n-dependent  portion  of  eq  (95)  may  be  used  in  eq  (92)  to  obtain 


sidered  in  this  section.  Curves  A,  B,  and  C  in  figure  7  show  the  temperature  cal- 
culated at  the  center  of  the  heat  source  as  a  function  of  the  number  of  terms  in 
the  sum,  N ,  for  the  top  of  the  first,  second,  and  third  layers,  respectively.  It  is 
clear  from  the  curve  in  figure  7a  that  the  calculation  of  the  surface  temperature  (at 
z  =  Qi)  may  require  up  to  at  least  350  terms  in  the  sum  while,  from  figures  7b  and 
7c,  the  temperature  below  the  surface  may  need  only  20  terms  to  be  retained  in  the 
sum.  From  the  curve  depicted  in  figure  7a,  it  may  be  argued  that  the  number  of 
terms  needed  to  adequately  represent  a  feature  size  of  Ax  (typical  of  the  lateral  size 
of  a  heat  source)  should  be  on  the  order  of  Lx/ Ax.  For  the  present  case,  this  ratio 


(98) 


where  the  dependence  of  T^(x,  y,  z)  on  N ,  the  upper  limit  of  summation,  is  con- 
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is  about  200.  Deeper  into  the  structure,  the  heat  flow  has  caused  the  approximate 
size  of  the  heat  current  to  spread  out  and  hence  fewer  terms  are  necessary.  This  is 
borne  out  by  the  behavior  of  the  temperature  in  the  figure. 
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GENERAL  DISCUSSION  OF  THE  TXYZ  PROGRAM 


The  annotated  listing  of  the  program  is  contained  in  the  appendix  of  this  report.  In 
the  present  section,  several  aspects  of  the  program  will  be  discussed.  The  purpose 
of  this  discussion  is  to  present  the  user  with  information  concerning  the  implemen- 
tation and/or  modification  of  the  program. 

The  first  item  is  the  format  in  which  input  data  are  to  be  entered.  As  listed 
in  the  program,  data  are  to  be  entered  in  a  fixed-field  format.  However,  many 
FORTRAN  compilers  support  free-field  or  list-directed  read  statements.  Free-field 
read  statements  are  convenient  as  they  allow  the  user  to  change  input  data  without 
having  to  be  concerned  with  the  tedious  process  of  lining  up  of  the  data  required 
by  fixed-field  read  statements.  Hence,  if  the  user's  FORTRAN  compiler  allows  for 
free-field  read  statements,  it  is  strongly  suggested  that  the  program  be  edited  such 
that  all  read  statements  are  of  this  form. 

The  next  item  is  concerned  with  the  situation  of  the  hyperbolic  functions  being 
or  not  being  built-in  functions.  If  these  functions  are  contained  in  the  computer's 
relocatable  (or  object  code)  mathematics  library,  then  there  are  no  changes  required 
in  the  code.  If,  however,  these  are  not,  then  it  is  necessary  to  write  separate  function 
subroutines  for  the  hyperbolic  functions.  As  the  hyperbolic  functions  are  easily 
generated  from  the  exponential  function,  this  is  a  straightforward  process. 

Next,  it  is  important  to  discuss  the  number  of  heat  sources  and  the  dimension 
statements  used  in  the  program.  In  its  present  form,  TXYZ  will  allow  up  to  twenty 
heat  sources.  This  was  deemed  sufficient  for  most  problems  of  interest.  However, 
should  the  user  want  to  use  more  than  this  number  of  heat  sources  (especially 
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for  the  case  of  a  number  of  nonuniform  heat  sources),  it  is  necessary  to  change 
the  corresponding  dimension  statements.  This  is  commented  upon  in  the  program 
and  the  change  is  relatively  straightforward.  The  only  other  comment  concerning 
dimension  statements  has  to  do  with  the  number  of  points  at  which  the  temperature 
is  to  be  calculated  and  the  maximum  number  of  terms  to  be  used  in  the  calculation 
of  the  temperature  using  eq  (92).  The  program  allows  up  to  501  as  the  maximum 
number  for  both  of  these  quantities.  This  was  considered  to  be  adequate  for  the 
resolution  of  the  spatial  variation  of  the  temperature  and  for  the  calculation  of  the 
temperature  (even  for  the  smallest  physically  realizable  heat  source  on  a  structure  of 
typical  lateral  dimensions).  Hence,  it  is  felt  that  501  is  a  good  maximum  value  which 
gives  sufficient  detail  without  excessive  use  of  computer  memory  space.  Should 
memory  allocation  be  a  problem,  it  may  become  necessary  to  reduce  these  numbers 
in  the  dimension  statements.  It  is  important  to  carefully  evaluate  any  reduction  in 
the  maximum  number  of  terms  in  the  Fourier  representation  of  the  temperature 
in  light  of  the  slow  convergence  of  the  surface  temperature  for  small  heat  sources. 
Finally,  in  regard  to  changing  the  dimension  statements,  it  is  important  to  change 
the  dimension  statements  in  not  only  the  main  program  but  also  in  the  function 
subroutines  as  the  dimensioned  variables  are  declared  in  common  and  hence  must 
have  the  same  dimensions  everywhere  (in  the  main  program  and  the  function 
subroutines)  to  avoid  any  register  difficulties. 

A  copy  of  this  FORTRAN  program  may  be  obtained  from  the  author  by  sending  a 
letter  of  request  and  a  computer  tape. 
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SELECTED  PROGRAM  EXAMPLES 


In  order  to  assist  the  user  in  the  implementation  of  the  TXYZ  program,  several 
examples  are  presented  in  this  section.  In  particular,  the  input  data  files  are  listed 
and  plots  of  the  temperature  profiles  are  presented  in  the  accompanying  figures. 
The  input  data  files  have  been  annotated  in  order  to  facilitate  familiarity  with  the 
reading  of  these  files.  This  is  especially  important  because  of  the  number  of  input 
variables  which  set  various  "switches"  inside  the  program.  The  corresponding  data 
files  will  be  sent  along  with  the  FORTRAN  listing  of  the  TXYZ  program  so  as  to 
simplify  the  implementation  of  the  program  and  its  subsequent  use. 

Two  examples  will  be  discussed  in  this  section.  The  first  is  for  the  situation  of  a 
single  small  heat  source  on  the  surface  of  the  top  layer.  The  second  example  is  that 
of  a  number  of  heat  sources  on  the  top  layer.  In  particular,  one  of  the  heat  sources 
is  nonuniform  and  its  construction  from  uniform  heat  sources  is  illustrated  in  the 
ccorresponding  input  data  file. 

The  first  example  is  that  of  a  1  mil  by  1  mil  uniform  heat  source  on  the  surface 
of  a  150  mil  by  200  mil  three-layer  structure.  The  specific  structure  is  that  of  15 
mils  of  sihcon  on  2  mils  of  die  attach  material  on  30  mils  of  a  substrate.  This  three- 
layer  structure  will  be  used  in  both  this  example  and  the  one  to  follow.  The  input 
data  file  for  this  first  example  is  listed,  with  annotation,  in  Table  1.  The  calculated 
temperature  is  presented  in  figure  8. 

The  second  example  is  for  the  same  three-layer  structure  as  used  in  the  first  example. 
However,  in  this  case,  there  are  several  heat  sources,  with  one  of  them  being 
nonuniform  and  built  up  from  uniform  heat  sources.  In  particular,  there  are  three 
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square  uniform  heat  sources  on  the  right  portion  of  the  surface  and  one  uniform 
heat  source  in  the  middle  of  the  surface.  The  nonuniform  heat  source  is  a  single 
rectangular  area.  The  surface  temperature  calculated  along  the  line  0  <  x  <  200, 
y  =  101  is  presented  in  figure  9.  The  input  data  file  is  presented  in  Table  2. 

In  addition  to  these  examples,  the  reader's  attention  is  drawn  to  the  recent  review 
[5]  where  this  program  is  used  in  the  thermal  evaluation  of  VLSI  packages  using 
test  chips. 

ACKNOWLEDGMENTS 

The  author  would  like  to  thank  Frank  F.  Oettinger  for  his  continued  interest  during 
the  course  of  this  work.  His  relentless  use  of  the  program  brought  many  of  the  early 
bugs  to  the  surface. 

The  author  would  also  like  to  thank  Stephen  E.  Ross  who  changed  the  code  into  a 
portable  FORTRAN  form.  His  assistance  in  getting  the  three-dimensional  plots  is 
also  appreciated. 


37 


REFERENCES 

1.  Kokkas,  A.  G.,  Thermal  Analysis  of  Multiple-Layer  Structures,  IEEE  Trans. 
Electron  Devices  ED-21,  674-681  (1974). 

2.  Carslaw,  H.  S.,  and  Jaeger,  J.  C,  Conduction  of  Heat  in  Solids,  Second  Edition 
(Oxford  University  Press,  1959). 

3.  Churchill,  R.  V.,  Fourier  Series  and  Boundary  Value  ProblemSj  Second  Edition 
(McGraw-Hill,  1963). 

4.  Ayres,  F.,  Jr.,  Matrices  (Schaum  Publishing  Co.,  1950). 

5.  Oettinger,  F.  F.,  Thermal  Evaluation  of  VLSI  Packages  Using  Test  Chips  -  A 
Critical  Review,  Solid  State  Technology  27  (2),  169-179  (1984). 


38 


I 


DISTRIBUTED 
HEATERS 


HEAT  SINK 


\    Ls  SEMICONDUCTOR 

W/ 

^  1  Lc   t           DIE  ATTACH 

Lj   t  SUBSTRATE 

Figure  1.  This  figure  presents  the  geometry  of  the  three-layer  structure  in  which 
the  steady-state  temperature  is  calculated. 
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Figure  2.  This  figure  shows  the  n-dependence  of  the  product  of  the  power  density 
function  and  the  cos  terms  in  eq  (92),  i.e.,  U{n)  cos{n7rx/Lx),  evaluated  at  the  center 
of  a  single  small  heat  source  which  is  in  the  form  of  a  thin  stripe  along  the  y- 
direction.  This  specific  form  of  the  product  is  used  as  there  is  complete  coverage 
along  the  y-direction  which  means  that  only  the  m  =  0  term  contributes.  The 
general  feature  seen  in  this  figure  is  that  this  function  does  not  asymptote  to  zero 
for  n  on  the  order  of  500.  This  general  kind  of  behavior  has  been  found  to  be  the 
case  for  all  situations  considered. 
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ri  (n.m.O) 


Figure  3.  The  top-layer  Fourier  coefl&cient,  ri(n,  m,  0),  is  shown  as  a  function  of  n 
and  m  and  at  =  0  for  a  three-layer  structure  where  all  three  layers  have  different 
thermal  conductivities.  The  specific  structure  is  for  15  mils  of  silicon  (/ci  =  0.00267) 
over  2  mils  of  die  attach  (/C2  =  0.00064)  over  30  mils  of  substrate  material  (/C3  = 
0.00999).  The  peak  value  of  ri  is  at  the  origin  of  the  n,m  plane. 
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Figure  4.  The  top-layer  Fourier  coefficient,  Ti{n,m,0),  is  shown  as  a  function  of  n 
and  m  and  at  2:  =  0  for  a  three-layer  structure  where  all  three  layers  have  the  same 
thermal  conductivity,  i.e.,  a  thick  one-layer  structure.  The  specific  example  is  for  47 
mils  of  silicon  (/c  =  0.00267).  As  with  the  previous  figure,  ri  is  peaked  at  the  origin 
of  the  n,  m  plane.  Both  this  and  the  previous  figure  are  meant  to  convey  qualitative 
features  of  the  Fourier  coefficient.  More  detailed  information  will  be  presented  in 
the  next  two  figures  where  specific  attention  will  be  focused  upon  just  the  n  axis 
(with  m  =  0)  behavior. 
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Figure  5.  This  jBgure  shows  the  dependence  of  the  three  Fourier  coefficients  on  n  at 
the  top  of  each  of  the  layers.  The  curve  denoted  by  A  is  Ti{n,m  =  0,z  =  0)  (for 
the  top  of  the  first  layer).  B  shows  the  Fourier  coefficient,  T2{n,m  =  0,z  =  —Li) 
(for  the  top  of  the  second  layer).  Finally,  the  curve  denoted  by  C  represents  the 
Fourier  coefficient,  Ts{n,m  =  0,z  =  —Li  —  L2),  i.e.,  at  the  top  of  the  third  layer. 
This  figure  and  the  next  represent  a  cross  section  along  the  n  axis  of  figures  of  the 
type  presented  in  figures  3  and  4. 
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Figure  6.  The  behavior  of  the  top  layer  Fourier  coefficient,  Ti{n,m  =  0,z  =  0), 
along  the  n  axis  is  presented  in  this  figure.  These  results  (curve  A  of  fig.  5)  are 
presented  in  a  log-log  plot  to  show  the  slow  convergence  z  =  0.  When  curves  B 
and  C  of  figure  5  are  plotted  in  the  same  manner,  they  fall  off  to  zero  much  more 
rapidly. 
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Figure  7.  This  figure  shows  the  calculated  temperature  as  a  function  of  the  number 
of  terms  used  in  the  sum  (eq  (98))  for  the  top  of  each  of  the  three  layers.  Curve  A  is 
Ti^(z  =  0)  while  curves  B  and  C  represent  T^{z  =  -Li)  and  T^{z  =  -Li  -  L2), 
respectively.  This  figure  clearly  shows  how  the  behavior  of  the  Fourier  coefficients 
for  each  of  the  layers  is  mirrored  in  the  corresponding  calculated  temperatures. 
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Figure  8.  This  figure  presents  the  calculated  temperature  along  the  midline  of  a  200 
mil  by  200  mil  rectangular  structure  with  a  single,  uniform  power  density  1  mil  by 
1  mil  heat  source  located  at  the  center  of  the  surface  of  the  top  layer.  Note  that 
only  the  temperature  in  the  region  of  the  heat  source  is  plotted,  i.e.,  from  90  mils 
to  110  mils.  Also,  note  that  the  temperature  falls  off  rapidly  away  from  the  heat 
source.  The  specific  structure  is  that  of  15  mils  of  silicon  on  2  mils  of  die  attach 
material  on  30  mils  of  substrate  material.  The  annotated  input  data  file  used  to 
obtain  this  temperature  data  is  to  be  found  in  Table  I. 
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Figure  9.  By  way  of  contrast  with  the  previous  figure,  the  calculated  temperature 
along  the  midline  of  a  200  mil  by  200  mil  three-layer  structure  is  presented  in  this 
figure.  However,  in  this  case,  two  nonuniform  heat  sources  and  one  uniform  heat 
source  are  used.  The  three-layer  structure  is  the  same  as  used  in  the  previous  figure. 
The  construction  of  the  nonuniform  heat  sources  from  a  number  of  uniform  heat 
sources  is  contained  in  the  annotated  input  data  file  contained  in  Table  11. 
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TABLE  I 

INPUT  DATA  FOR  SINGLE  HEAT  SOURCE 


INPUT  DATA 


DESCRIPTION  OF  INPUT  DATA 


150 
0.00267 
0.00064 
0.00999 
500 

90  0.1 


200 
15 
2 

30 
500 
1 

201 
0 

75.5 

0 

0 

1  0.093423 
99      1  75 

0 
0 

0  0 
0 


X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 
UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 
lEDGEX   (=0  FOR  SINGLE  POINT,  =1   FOR  NUMBER  OF  POINTS) 
NUMBER  OF  POINTS  ALONG  X,   FIRST  POINT,  STEP  INCREMENT 
lEDGEY  (=0  FOR  SINGLE  POINT,  =1   FOR  NUMBER  OF  POINTS) 
Y  POINT  FOR  CALCULATION 

lEDGEZ  (=0  FOR  SINGLE  POINT,  =1   FOR  NUMBER  OF  POINTS) 
Z  POINT  FOR  CALCULATION 

NUMBER  OF  HEAT  SOURCES  AND  POWER  DENSITY 

X  COORDINATE,  LENGTH  ALONG  X  AXIS,  Y  COORDINATE,  LENGTH 

ALONG  Y  AXIS 
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TABLE  II 

INPUT  DATA  FOR  MULTIPLE,  NONUNIFORM  HEAT  SOURCES 


INPUT  DATA 


DESCRIPTION  OF  INPUT  DATA 


200  200 
15  0.00267 
2  0.00064 
30  0.00999 
500  500 
1 

200      0  1 
0 

101 

0 

0 

7  0.001 

10      65      10  180 


30 
50 
100 
160 


45 
25 
2 


10  180 
10  180 
100  2 


20      10  20 


160      20      90  20 

160      20      170  20 

0 
0 

0  0 
0 


X  AND  Y  DIMENSIONS  OF  RECTANGULAR  STRUCTURE 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  TOP  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER 
THICKNESS  AND  THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER 
UPPER  SUMMATION  LIMITS  FOR  N  AND  M  SUMMATIONS 
lEDGEX   (=0  FOR  SINGLE  POINT,  =1   FOR  NUMBER  OF  POINTS) 
NUMBER  OF  POINTS  ALONG  X,  FIRST  POINT,  STEP  INCREMENT 
lEDGEY  (=0  FOR  SINGLE  POINT,  =1   FOR  NUMBER  OF  POINTS) 
Y  POINT  FOR  CALCULATION 

lEDGEZ  (=0  FOR  SINGLE  POINT,  =1  FOR  NUMBER  OF  POINTS) 
Z  POINT  FOR  CALCULATION 

NUMBER  OF  HEAT  SOURCES  AND  POWER  DENSITY 


X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #1 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #2 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #3 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #4 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  US 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #6 

X  COORDINATE,  LENGTH  ALONG  X  AXIS, 

ALONG  Y  AXIS  FOR  HEATER  #7 


COORDINATE,  LENGTH 

COORDINATE,  LENGTH 

COORDINATE,  LENGTH 

COORDINATE,  LENGTH 

COORDINATE,  LENGTH 

COORDINATE,  LENGTH 

COORDINATE,  LENGTH 
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MAIN  PROGRAM  LISTING 


C  ********************************************************** 

C  INSTRUCTIONS  FOR  INPUT  DATA  FILE  SETUP 

C  DATA  TO  BE  READ  BY  THIS  PROGRAM  ARE  SET  UP  AS  FOLLOWS: 
C 

C  RLX   (X  DIMENSION  OF  3  LAYER  STRUCTURE) 

C  RLY   (Y  DIMENSION  OF  3  LAYER  STRUCTURE) 

C  RLS   (THICKNESS  OF  TOP  LAYER) 

C  RKS   (THERMAL  CONDUCTIVITY  OF  TOP  LAYER) 

C  RLC   (THICKNESS  OF  MIDDLE  LAYER) 

C  RKC   (THERMAL  CONDUCTIVITY  OF  MIDDLE  LAYER) 

C  RLI   (THICKNESS  OF  BOTTOM  LAYER) 

C  RKI   (THERMAL  CONDUCTIVITY  OF  BOTTOM  LAYER) 

C  NUP  (UPPER  LIMIT  OF  N  SUM,  X  DIRECTION) 

C  MUP  (UPPER  LIMIT  OF  M  SUM,  Y  DIRECTION) 

C 

C  lEDGEX   (=0  FOR  DEFAULT,  =1  TO  ALTER  RANGE  OF  X  VALUES) 

C  IF  IEDGEX=1   THEN  THE  FOLLOWING  SEGMENT  OF  INPUT  IS  READ 

C  ILX  (THE  NUMBER  OF  INCREMENTS  IN  X  TO  BE  USED) 

C  (ILX+1=THE  NUMBER  OF  X  VALUES  TO  BE  USED) 

C  XI    (THE  VALUE  OF  THE  FIRST  POINT  IN  X) 

C  STEPX   (THE  INCREMENT  IN  X) 

C  ELSE  IEDGEX=0     READ  XI    (THE  CONSTANT  X  VALUE) 

C 

C  lEFGEY  (=0  FOR  DEFAULT,  =1  TO  ALTER  RANGE  OF  Y  VALUES) 

C  IF  IEDGEY=1   THEN  THE  FOLLOWING  SEGMENT  OF  INPUT  IS  READ 

C  ILY  (THE  NUMBER  OF  INCREMENTS  IN  Y  TO  BE  USED) 

C  (ILY+1=THE  NUMBER  OF  Y  VALUES  TO  BE  USED) 

C  Y1   (THE  VALUE  OF  THE  FIRST  POINT  IN  Y) 

C  STEPY   (THE  INCREMENT  IN  Y) 

C  ELSE  IEDGEX=0     READ  XI    (THE  CONSTANT  X  VALUE) 

C 

C  lEDGEZ  (=0  FOR  DEFAULT,  =1  TO  ALTER  RANGE  OF  Z  VALUES) 

C  IF  IEDGEZ=1   THEN  THE  FOLLOWING  SEGMENT  OF  INPUT  IS  READ 

C  ILZ  (THEN  NUMBER  OF  INCREMENTS  IN  Z  TO  BE  USED) 

C  (ILZ+1=THE  NUMBER  OF  Z  VALUES  TO  BE  USED) 

C  Z1    (THE  VALUE  OF  THE  FIRST  POINT  IN  Z) 

C  STEPZ  (THE  INCREMENT  IN  Z) 

C  ELSE  IEDGEX=0     READ  XI    (THE  CONSTANT  X  VALUE) 

C 

C  NOTE:  AS  THE  CALCULATION  TAKES  THE  Z  VARIABLE  TO  BE  ZERO  OR 
C  NEGATIVE,  ENTER  Z1   AND  STEPZ  AS  POSITIVE  QUANTITIES. 

C  THE  PROGRAM  CHANGES  THE  SIGN  OF  Z1   AND  STEPZ  TO  MAKE 

C  THE  CALCULATION  FOR  NEGATIVE  Z. 

C 

C  NSOUR  (NUMBER  OF  HEAT  SOURCES,  UP  TO  20) 

C  PO   (POWER  DENSITY) 

C  XSOUR(I)    (X  COORDINATE  OF  ORIGIN  OF  1ST  SOURCE) 
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C  YSOUR(I)    (Y  COORDINATE  OF  ORIGIN  OF  1ST  SOURCE) 

C  RLXSOUR(I)    (LENGTH  ALONG  X  DIRECTION  OF  1ST  SOURCE) 

C  RLYSOUR(I)    (LENGTH  ALONG  Y  DIRECTION  OF  1ST  SOURCE) 

C  REMAINING  HEAT  SOURCES  WITH  SAME  INPUT  STRUCTURE  AS 
C 

C  ***************************************************** 

c 

DIMENSION  X(501),Y(501),Z(501),COSYT(501) 
DIMENSION  ARUZER(501,501 ),ARFUNZ(501,501 ) 
DIMENSION  XSOUR(20),YSOUR(20),RLXSOR(20),RLYSOR(20) 
COMMON  RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI 
COMMON  NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR 

C 

c  ***************************************************** 

c 

C  DOCUMENTATION  AND  BACKGROUND 

C  THIS  PROGRAM  CALCULATES  THE  SURFACE  TEMPERATURE  T(X,Y,Z)   DUE  TO  DC 

C  POWER  INPUTS  ONLY.     THE  TEMPERATURE  IS  THE  TEMPERATURE  RELATIVE 

C  TO  THE  AMBIENT.     THE  SPECIFIC  EQUATIONS  USED  ARE  GIVEN  IN 

C  EQUATIONS   (13)-(23),  WITH  S=0   (STEADY-STATE  CONDITION), 

C  IN  THE  PAPER  BY  KOKKAS   (REF:  "THERMAL  ANALYSIS  OF  MULTIPLE- 

C  LAYERED  STRUCTURES"  BY  ACHILLES  G.  KOKKAS,  IEEE  TRANS.   ELEC.  DEV. 

C  VOL.  ED-21,  NO.  11,  674-681  (1974)). 

C  VARIABLES  USED:   THE  VARIABLES  LISTED  AS  REAL  IN  THE  ABOVE 

C  ARE  THE  FOLLOWING-KS,KC,  AND  KI  ARE  THE  THERMAL  CONDUCTIVITIES 

C  OF  THE  SEMICONDUCTOR,  CONDUCTOR,  AND  INSULATOR,  RESPECTIVELY. 

C  LX  AND  LY  ARE  THE  LATERAL  DIMENSION  OF  THE  CHIP  WHILE 

C  LS,LC,  AND  LI  ARE  THE  THICKNESSES  OF  THE  SEMICONDUCTOR,  THE 

C  CONDUCTOR,  AND  THE  INSULATOR,  RESPECTIVELY. 

C  IMPORTANT  NOTE:     WHILE  THE  VARIABLES  HAVE  THE  NOTATION  WHICH 

C  SEEMS  TO  IMPLY  A  SEMICONDUCTOR,  A  CONDUCTOR,  AND  AN  INSULATOR, 

C  THESE  REFER  TO  THE  WAY  IN  WHICH  KOKKAS  FORMULATED  THE  PROBLEM. 

C  THE  THING  TO  KEEP  IN  MIND  IS  THAT  THE  TOP  LAYER  HAS  (LS,KS), 

C  THE  MIDDLE  LAYER  HAS   (LC,KC),  AND  THE  BOTTOM  LAYER  HAS   (LI,KI)  . 

C  IT  IS  NOT  NECESSARY  THAT  THEY  BE  WHAT  THEY  SEEM  TO  BE,  THEY  ARE 

C  DETERMINED  BY  THE  RESPECTIVE  THICKNESSES,  L,  AND  CONDUCTIVITIES, K . 
C 

C  ************************************************************* 

c 

1  FORMAT(IHI) 

2  FORMATdX, 'STEADY-STATE  THERMAL  ANALYSIS  CALCULATION  USING  EQS . 
1    (13)-(23)   OF  KOKKAS'/) 

3  FORMATdX, 'THERMAL  CONDUCTIVITIES  AND  LAYER  THICKNESSES') 

4  FORMATdX, 'KS=   ',F10.8,'   KC=   ',F10.8,'   KI=  ',F10.8) 

5  FORMATdX, 'LS=   ',F10.5,'   LC=   ',F10.5,'   LI=  ',F10.5) 

6  F0RMAT(/1X, 'UPPER  SUMMATION  LIMITS   ',2X,'  NUP=',I5, 
(  1    '  MUP=',I5/) 

7  F0RMAT(//1 X, 'NUMBER  OF  HEAT  SOURCES= ' ,15) 

8  F0RMAT(/1X, 'COORDINATES,  LENGTHS,  AND  WIDTHS  OF  HEAT  SOURCES'/) 

9  FORMATdX, 'HEAT  SOURCE  '  ,7X,  '  XSOUR '  ,1  OX,  '  YSOUR '  ,9X,  '  LXSOUR '  , 
1  9X,'LYS0UR'/) 

10  F0RMAT(7X,I3,5X,F10.5,5X,F10.5,5X,F10.5,5X,F10.5) 
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11  FORMATdX, 'POWER  DENSITY=' ,F1  2.6) 

12  F0RMAT(/1X, 'CALCULATING   ',13,'  X  POINTS  WITH  A  FIRST  POINT  OF  ' 
1     ,F5.1,'   AND  A  STEP  SIZE  OF  ',F5.1) 

13  F0RMAT(/1X,'THE  CONSTANT  X  COORDINATE  IS  ',F5.1) 

14  F0RMAT(/1X, 'CALCULATING   ',13,'   Y  POINTS  WITH  A  FIRST  POINT  OF  ' 
1     ,F5.1,'   AND  A  STEP  SIZE  OF  ',F5.1) 

15  F0RMAT(/1X,'THE  CONSTANT  Y  COORDINATE  IS  ',F5.1) 

16  F0RMAT(/1X, 'CALCULATING  ',13,'   Z  POINTS  WITH  A  FIRST  POINT  OF  ' 
1     ,F5.1,'   AND  A  STEP  SIZE  OF  ',F5.1) 

17  F0RMAT(/1X,'THE  CONSTANT  Z  COORDINATE  IS  ',F5.1) 
27  F0RMAT(1X,'LX=  ',F7.2,3X,'  LY=  ',F7.2) 

22  F0RMAT(1X,3F8.2,F10.4) 

31  F0RMAT(1X,6I7) 

51  F0RMAT(I4,3X,I4) 

52  F0RMAT(F10.5,3X,F10.5) 

53  FORMAT(II) 

54  FORMAT(I4,3X,F10.5,3X,F10.5) 

55  FORMATCFIO.S) 

56  F0RMAT(I2,3X,F10.5) 

57  F0RMAT(F10.5,3X,F10.5,3X,F10.5,3X,F10.5) 

C 

C  ********************************************** *********^ 

C  DATA  INPUT  SECTION   (DATA  READ  FROM  F0R010) 

C  ****************************************************************** 

c 

IXFLG  =  0 
lYFLG  =  0 
IZFLG  =  0 
ILX=0 
ILY=0 
ILZ=0 
STEPX=1.0 
STEPY=1.0 
STEPZ=-1 .0 
WRITE(11,1 ) 
WRITE(11,2) 
WRITE(11,3) 

C  RLX  AND  RLY  ARE  THE  X  AND  Y  DIMENSIONS  OF  THE  RECTANGULAR  STRUCTURE 

READ(10,52)RLX,RLY 
WRITE(11,27)RLX,RLY 

C  RLS  IS  THE  THICKNESS  AND  RKS  IS  THE  THERMAL  CONDUCTIVITY  OF  THE  TOP  LAYER 

READ(10,52)RLS,RKS 

C  RLC  IS  THE  THICKNESS  AND  RKC  IS  THE  THERMAL  CONDUCTIVITY  THE  MIDDLE  LAYER 

READ(10,52)RLC,RKC 

C  RLI  IS  THE  THICKNESS  AND  RKI  IS  THE  THERMAL  CONDUCTIVITY  THE  BOTTOM  LAYER 

READ(10,52)RLI,RKI 

WRITE(11,5)RLS,RLC,RLI 

WRITE (11, 4) RKS, RKC, RKI 
C  NUP  IS  THE  UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE  INDEX  N  (X-DIR) 

C  MUP  IS  THE  UPPER  LIMIT  OF  THE  SUMMATION  OVER  THE  INDEX  M  (Y-DIR) 

READ(10,51)NUP,MUP 

WRITE(11,6)NUP,MUP 
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NUP  =  NUP  +  1 
MUP  =  MUP  +  1 

IF  (NUP. GT. 501 .OR. MUP. GT. 501)  GO  TO  3999 
888     FORMATCIX,' YOUR  UPPER  LIMIT  OF  SUMMATION  IS  TOO  LARGE.  TRY  AGAIN') 
READ(10,53)IEDGEX 
C  SET  IEDGEX=0  TO  GET  CONSTANT  POINT  FOR  X 

C  SET  IEDGEX=1   TO  SET  THE  LIMITS  OF  THE  CALCULATION  ALONG  X 

IF(IEDGEX.EQ.O)  GO  TO  115 

READ (10,54) ILX,X1,STEPX 

WRITE(11,12)ILX,X1,STEPX 

GO  TO  125 
115  READ(10,55)X1 

WRITE(11,13)X1 

IXFLG  =  1 
125  READ(10,53)IEDGEY 
C  SET  IEDGEY=0  TO  GET  CONSTANT  POINT  FOR  Y 

C  SET  IEDGEY=1   TO  SET  THE  LIMITS  OF  THE  CALCULATION  ALONG  Y 

IF(IEDGEY.EQ.O)  GO  TO  135 

READ  (10,54) ILY,Y1,STEPY 

WRITE (11,1 4) ILY,Y1,STEPY 

GO  TO  150 
135  READ(10,55)Y1 

WRITE(11,15)Y1 

lYFLG  =  1 
150  READ(10,53)IEDGEZ 
C  SET  IEDGEZ=0  TO  GET  CONSTANT  POINT  FOR  Z 

C  SET  IEDGEZ=1   TO  SET  THE  LIMITS  OF  THE  CALCULATION  ALONG  Z 

IF   (lEDGEZ.EQ.O)  GO  TO  160 

READ  (10,54)ILZ,Z1,STEPZ 

WRITEd  1,16)ILZ,Z1  ,STEPZ 

Z1  =-1.0*Zl 

STEPZ=-1 .0*STEPZ 

GOTO  175 
160  READ(10,55)Z1 

WRITE(11,17)Z1 

Z1  =  -1.0*Z1 

IZFLG  =  1 

ILX  =  ILX  +  1 

ILY  =  ILY  +  1  . 

ILZ  =  ILZ  +  1 
175  READ(10,56)NS0UR,P0 
C  PO  IS  THE  POWER  DENSITY,  ASSUMED  UNIFORM  FOR  ALL  HEATERS 

C  NSOUR  IS  THE  TOTAL  NUMBER  OF  HEATING  ELEMENTS  ON  THE  SURFACE  OF  THE 

C  THE  TOP  LAYER  (UP  TO  20  ARE  POSSIBLE) 

WRITE(11,7)NS0UR 

WRITE(11,11)  PO 

WRITE(11,8) 

WRITE(11,9) 

C  THE  NEXT  LOOP  READS  IN  THE  COORDINATES  OF  THE  ORIGIN  OF  THE 

C  HEATING  ELEMENTS  ALONG  WITH  THEIR  LENGTHS  AND  WIDTHS 

DO  100  1=1, NSOUR 

READ(10,57)XS0UR(I),RLXS0R(I),YS0UR(I),RLYS0R(I) 
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C  XSOUR(I)   IS  THE  X  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 

C  RLXSOR(I)   IS  THE  LENGTH  OF  THE  I-TH  HEATER  ALONG  THE  X  DIRECTION 

C  YSOUR(I)   IS  THE  Y  COORDINATE  OF  THE  ORIGIN  OF  I-TH  HEATER  ELEMENT 

C  RLYSOR(I)   IS  THE  LENGTH  OF  THE  I-TH  HEATER  ALONG  THE  Y  DIRECTION 


WRITE(11,10)I,XS0UR(I),YS0UR(I),RLXS0R(I),RLYS0R(I) 
100  CONTINUE 


WRITE(11,1  ) 

C 

C  ******** ***************************************** 
C 

C  END  OF  DATA  INPUT  SECTION 

C  END  OF  INPUT  SECTION.     THE  THERMAL  CONDUCTIVITIES  OF  THE 

C  SEMICONDUCTOR,  CONDUCTOR,  AND  INSULATOR,  THE  THICKNESSES  OF  THE 

C  SEMICONDUCTOR,  CONDUCTOR,  AND  INSULATOR,  THE  X  AND  Y  DIMENSIONS 

C  OF  THE  CHIP,  THE  NUMBER  OF  HEATING  SOURCES  AND  THEIR  X,Y  AND 

C  LENGTH  AND  WIDTH  HAVE  BEEN  ENTERED. 

C 

C  ***************************************************************** 

c 


180  PI=3. 14159265 

P04LK  =  4.0  *  PO  /  (RLX*RLY*RKS) 
PILX  =  PI  /  RLX 
PILY  =  PI  /  RLY 

C 

C  ****************************************************************** 

C  CALCULATE  THE  X(0:INTX),  Y(0:INTY),  AND  Z(0:INTZ)  ARRAYS 

C  ****************************************************************** 

c 

X(1)=X1 

DO  200  1=2, ILX 
X(I)=X(I-1)+STEPX 
200  CONTINUE 
Y(1)=Y1 

DO  220  1=2, ILY 
Y(I)=Y(I-1)+STEPY 
220  CONTINUE 
Z(1)=Z1 

DO  240  I=2,ILZ 
Z(I)=Z(I-1)+STEPZ 
240  CONTINUE 


C 

C  *************************************************************** 
C 

C  BEGIN  CALCULATION  OF  T(X,Y,Z) 

C  THE  SUBROUTINES  USED  IN  THE  CALCULATION  ARE: 

C  1)  UZERO(N,M)  -  CALCULATES  THE  FOURIER  COSINE  TRANSFORM  OF  THE 
C  FUNCTION,  U(X,Y),  THE  POWER  DENSITY  FUNCTION  FOR  ALL  OF  THE 

C  HEAT  SOURCES. 

C  2)   FUNZ(N,M,Z)  -  CALCULATES  THE  Z-DEPENDENT  PORTION  OF  THE  SUM 
C  REMEMBERING  THAT  THIS  IS  A  FUNCTION  OF  THE  SUMMATION 

C  INDICES  (N,M). 

C 
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C  ******************************************************* 

c 
c 

c  ***************************************************************** 

C  CALCULATE  THE  FOURIER  COMPONENTS  OF  THE  HEAT  SOURCES,  U(N,M) 

C  ***************************************************************** 

c 

DO  300  MM=1,MUP 
M  =  MM  -  1 
DO  250  NN=1,NUP 
N  =  NN  -  1 

ARUZER(NN,MM)=UZERO(N,M) 
250  CONTINUE 
300  CONTINUE 

C 

c  ***************************************************************** 

c 

C  END  OF  U(N,M)   CALCULATION  AND  BEGINNING  OF  MAJOR  LOOP  FOR  Z 
C 

c  ***************************************************************** 

c 

DO  3000  IZ=1,ILZ 

c 

c  ***************************************************************** 

C  CALCULATE  THE  Z  DEPENDENT  POTION,  I.E.,  FUNZ(N,M,Z) 

C  ***************************************************************** 

c 

DO  400  MM=1,MUP 
M  =  MM  -  1 
DO  350  NN=1,NUP 
N  =  NN  -  1 

ARFUNZ(NN,MM)=FUNZ(N,M,Z(IZ))*ARUZER(NN,MM) 
350  CONTINUE 
400  CONTINUE 

C 

DO  3000  IY=1,ILY 
DO  700  MM=1,MUP 
M  =  MM  -  1 

COSYT(MM)=COS(FLOAT(M)*Y(IY)*PILY) 
700  CONTINUE 

C 

DO  3000  IX=1,ILX 
SUM=0.0 

DO  1900  MM=1,MUP 

M  =  MM  -  1 

DO  1700  NN=1,NUP 

N  =  NN  -  1 

NDN=0 

NDM=0 

IF  (N.EQ.O)  NDN=1 
IF   (M.EQ.O)  NDM=1 

TOP  =  ARFUNZ(NN,MM)   *  COS ( FLOAT(N) *X (IX) *P ILX)   *  COSYT(MM) 
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B0TT0M=(NDN+1 )*(NDM+1) 

TSUM=TOP/BOTTOM 

SUM=SUM+TSUM 
1700  CONTINUE 
1900  CONTINUE 

TEMP  =  P04LK  *  SUM 

WRITE  (20,22)X(IX),Y(IY),Z(IZ),TEMP 
3000  CONTINUE 

GO  TO  4000 

3999  WRITE(6,888) 

4000  STOP 
END 

C  **************************************************************^ 

C  END  OF  THE  MAIN  PROGRAM 

C  ******************************************************* 
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FUNCTION  UZERO  LISTING 


FUNCTION  UZERO(N,M) 

DIMENSION  XSOUR(20),YSOUR(20),RLXSOR(20),RLYSOR(20) 
COMMON  RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI 
COMMON  NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR 

C 

C  ******************* ************************************* 

C  DESCRIPTION  OF  THE  FUNCTION  UZERO(N,M) 

C  THIS  FUNCTION  CALCULATES  THE  DOUBLE  FOURIER  COSINE  TRANSFORM 

C  OF  THE  POWER  DENSITY  FUNCTION,  U(X,Y).     THIS  IS  THE  TRANSFORM 

C  FOR  ALL  OF  THE  HEAT  SOURCES.     THE  ASSUMPTION  IS  MADE  THAT  THE 

C  POWER  DENSITY  IS  UNIFORM  AND  EQUAL  TO  UNITY  OVER  THE  SURFACE 

C  OF  THE  HEATING  ELEMENTS.     THAT  IS, 

C  U(X,Y)=1     (XSOUR(I)<=X<=XSOUR(I)+LXSOUR(I)  AND 

C  YSOUR(I)<=Y<=YXOUR(I)+LYSOUR(I)  ). 

C  U(X,Y)=0  OTHERWISE. 

C  UNDER  THESE  CONDITIONS,  IT  IS  POSSIBLE  TO  ANALYTICALLY  EVALUATE 

C  THE  DOUBLE  INTEGRAL  FOR  EACH  HEATING  ELEMENT.     AS  THE  HEATING 

C  ELEMENTS  ARE  ASSUMED  TO  BE  INDEPENDENT,  THE  CONTRIBUTION  FROM 

C  EACH  ELEMENT  MAY  BE  ADDED  TO  OBTAIN  THE  U(N,M)    FOR  ALL. 

C 

C  *************************************************************** 
C 

PI=3.U159265 

UZERO=0.0 

DO  500  I=1,NS0UR 

IF(N.EQ.O)   GO  TO  100 

TERMX  =  SIN(FLOAT(N)*PI*(XSOUR(I)  +  RLXSOR (I) ) /RLX) 
1  -  SIN(FLOAT(N)*PI*XSOUR(I)/RLX) 

TERMX=TERMX*RLX/ (FLOAT(N)*PI) 
GO  TO  150 
100  TERMX=RLXSOR(I) 
150  IF(M.EQ.O)   GO  TO  200 

TERMY  =  SIN(FLOAT(M)*PI*(YSOUR(I)  +  RLYSOR (I) ) /RLY) 
1  -  SIN(FLOAT(M)*PI*YSOUR(I)/RLY) 

TERMY=TERMY*RLY/(FLOAT(M)*PI) 
GO  TO  250 
200  TERMY=RLYSOR(I) 
250  TERMI=TERMX*TERMY 
UZERO=UZERO+TERMI 
500  CONTINUE 
RETURN 
END 

C  ***************************************************************** 
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FUNCTION  FUNZ  LISTING 


FUNCTION  FUNZ(N,M,Z) 

DIMENSION  XS0UR(20),YS0UR(20),RLXS0R(20),RLYS0R(20) 
COMMON  RKS,RKC,RKI,RLX,RLY,RLS,RLC,RLI 
COMMON  NSOUR,XSOUR,YSOUR,RLXSOR,RLYSOR 

C 

C  ******************************************* 
C 

C  DESCRIPTION  OF  THE  FUNCTION  FUNZ(N,M,Z) 

C  THIS   FUNCTION  IS  USED  TO  CALCULATE  THE  Z  DEPENDENT  PART  OF  THE 

C  FUNCTION  USING  THE  S=0  VERSIONS  OF  EQUATIONS  (15)-(17)  OF  KOKKAS. 

C  THESE  ARE  USED  IN  CONJUNCTION  WITH  EQUATIONS   (18)-(22)  ALSO  FOR 

C  S=0  (STEADY-STATE  CONDITION).     THE  SPECIFIC  FORM  OF  FUNZ  IS 

C  DETERMINED  BY  THE  VALUE  OF  1,  I.E.,  IF  Z  FALLS  IN  THE  TOP,  MIDDLE 

C  OR  BOTTOM  LAYER  OF  THE  STRUCTURE.     IN  ADDITION,  SPECIAL  CARE  IS 

C  TAKEN  AS  TO  EVALUATE  FUNZ  FOR  THE  CASE  WHERE  GAMMA=0  AS  THESE  ARE 

C  SIMPLE,  BUT  THE  COMPUTER  DOES  NOT  KNOW  HOW  TO  EVALUATE  LIMITS. 

C  IN  ADDITION,  THE  FORM  OF  THE  THREE  SOLUTIONS  HAS  BEEN  CHANGED  TO 

C  GET  RID  OF  THE  ARTIFICIAL  OVERFLOW  PROBLEMS  COMING  FROM  THE 

C  HYPERBOLIC  FUNCTIONS,  COSH  AND  SINH,  FOR  THE  CASES  WHERE  THE 

C  ARGUMENTS  BECOME  LARGE. 

C 

c  ***************************************************** 

c 

PI=3. 14159265 

GAMMA=SQRT((FL0AT(N)*PI/RLX)**2  +  ( FLOAT(M) *PI/RLY) **2) 

VS=GAMMA*RLS 

VC=GAMMA*RLC 

VI=GAMMA*RLI 

VT=GAMMA*(RLS+Z) 

VM=GAMMA*(RLS+RLC+Z) 

VB=GAMMA*(RLS+RLC+RLI+Z) 

B0T1=TANH(VS)*TANH(VI) 

B0T1=B0T1+(RKI/RKC)*TANH(VS)*TANH(VC) 

B0T2=(RKC/RKS)*TANH(VI)*TANH(VC)+(RKI/RKS) 

GFUNC=1 .0/(B0T1+B0T2) 

AZ=ABS(Z) 

IF  (AZ.GT.RLS)  GO  TO  500 
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C  ************************************************************* 
C  TOP  LAYER  CALCULATION 

C  THIS  PORTION  IS  THE  TOP  LAYER  CALCULATION  WHICH  IS  DEFAULTED 

C  TO  IF  Z  FALLS  INTO  THE  TOP  LAYER 

C 

IF  (GAMMA. EQ. 0.0)   GO  TO  100 
TERMS1=TANH(VI)+(RKI/RKC)*TANH(VC) 

TERMS2=(RKC/RKS)*TANH(VT)*(TANH(VI)*TANH(VC)+(RKI/RKC) ) 

TERMS=TERMS1+TERMS2 

IF  (Z.EQ.0.0)  GO  TO  90 

IF   (VS.GT.5.0.AND.VT.GT.5.0)  GO  TO  80 

IF  (VS.LT.5.0)  GO  TO  10 

C1=2.0*EXP(-VS) 

GO  TO  20 
10      C1=1 .0/COSH(VS) 
20  CONTINUE 

IF  (VT.LT.5.0)  GO  TO  30 

C2=0.5*EXP(VT) 

GO  TO  40 
30  C2=C0SH(VT) 
40  CONTINUE 

FUNZ=GFUNC*TERMS*C1*C2/GAMMA 

RETURN 

80       FUNZ=GFUNC*TERMS*EXP(GAMMA*Z) /GAMMA 
RETURN 

90  FUNZ=GFUNC*TERMS/GAMMA 
RETURN 

100  FUNZ=(RLS+Z)+(RKS/RKC)*RLC+(RKS/RKI)*RLI 

RETURN 
500  TOTAL=RLS+RLC 

IF  (AZ.GT. TOTAL)  GO  TO  1500 
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C  ********************************************************** 

C  MIDDLE  LAYER  CALCULATION 

C  THIS  IS  THE  MIDDLE  LAYER  CALCULATION  WHICH  IS  DEFAULTED  TO  IF 

C  Z  FALLS  INTO  THIS  DOMAIN  OF  DEPTHS 

C 

IF  (GAMMA. EQ. 0.0)   GO  TO  1000 
TERMC=TANH(VI)+(RKI/RKC)*TANH(VM) 

IF  (VS.GT.5.0.AND.VC.GT.5.0.AND.VM.GT.5.0)  60  TO  800 

IF  (VS.LT.5.0)  GO  TO  250 

C1=2.0*EXP(-VS) 

60  TO  260 
250  C1=1.0/C0SH(VS) 
260  CONTINUE 

IF  (VC.LT.5.0)  GO  TO  270 

C2=2.0*EXP(-VC) 

GO  TO  280 
270     C2=1 .0/COSH(VC) 
280  CONTINUE 

IF  (VM.LT.5.0)  GO  TO  320 

C3=0.5*EXP(VM) 

GO  TO  330 
320  C3=C0SH(VM) 
330  CONTINUE 

FUNZ=GFUNC*TERMC*C1*C2*C3/GAMMA 

RETURN 

800     FUNZ=6FUNC*TERMC*2.0*EXP(6AMMA*Z) /GAMMA 
RETURN 

1000  FUNZ=(RKS/RKI)*RLI+(RKS/RKC)*(RLS+RLC+Z) 
RETURN 
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************ ****************************************** 

BOTTOM  LAYER  CALCULATION 
THIS  IS  THE  BOTTOM  LAYER  CALCULATION  WHICH  IS  USED  IF  Z  FALLS 
INTO  THE  BOTTOM  LAYER 


1500  IF  (GAMMA. EQ. 0.0)  GO  TO  2000 

IF  (VS.GT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0.AND. 
1  VB.GT.5.0)  GO  TO  1900 

IF(VB.GT.5.0.AND.VS.GT.5.0. AND.VC.LT.5.0. AND.VI.LT.5.0)G0  TO  2100 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2200 

IF(VB.GT.5.0.AND.VS.LT.5.0. AND.VC.LT.5.0. AND. VI. GT. 5. 0)G0  TO  2300 

IF(VB.GT.5.0.AND.VS.GT.5.0.AND.VC.GT.5.0.AND.VI.LT.5.0)GO  TO  2400 

IF(VB.GT.5.0.AND.VS.GT.5.0. AND.VC.LT.5.0. AND. VI. GT. 5. 0)G0  TO  2500 

IF(VB.GT.5.0.AND.VS.LT.5.0.AND.VC.GT.5.0.AND.VI.GT.5.0)GO  TO  2600 

IF  (VS.LT.5.0)  GO  TO  1550 

C1=2.0*EXP(-VS) 

GO  TO  1560 
1550  C1=1 .0/COSH(VS) 
1560  CONTINUE 

IF  (VC.LT.5.0)  GO  TO  1570 

C2=2.0*EXP(-VC) 

GO  TO  1580 
1570  C2=1 .0/COSH(VC) 
1580  CONTINUE 

IF  (VI.LT.5.0)  GO  TO  1590 

C3=2.0*EXP(-VI) 

GO  TO  1600 
1590  C3=1 .0/COSH(VI) 
1600  CONTINUE 

C4=SINH(VB) 

FUNZ=GFUNC*C1*C2*C3*C4/GAMMA 
RETURN 

1900  FUNZ=GFUNC*4.0*EXP(GAMMA*Z) /GAMMA 
RETURN 

2000  FUNZ=(RKS/RKI)*(RLS+RLC+RLI+Z) 
RETURN 

2100  FUNZ=GFUNC*EXP(VB-VS) /(GAMMA*COSH(VC)*COSH(VI)) 
RETURN 

2200  FUNZ=GFUNC*EXP(VB-VC)/(GAMMA*COSH(VS)*COSH(VI)) 
RETURN 

2300  FUNZ=GFUNC*EXP(VB-VI) /(GAMMA*COSH(VS)*COSH(VC) ) 
RETURN 

2400  FUNZ=GFUNC*2.0*EXP(VB-VS-VC) /(GAMMA*COSH(VI) ) 
RETURN 

2500  FUNZ=GFUNC*2.0*EXP(VB-VS-VI) /(6AMMA*C0SH(VC) ) 
RETURN 

2600  FUNZ=GFUNC*2.0*EXP(VB-VC-VI) /(GAMMA*C0SH(VS)) 
RETURN 
END 
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